Fishery resource management with migratory prey harvesting in two zones- delay and stochastic approach

In this work, we looked at a two-zone aquatic habitat where both prey and predators can access the zones. The prey alternates between two zones at random. The growth of prey in the absence of a predator is believed to be logistic in each zone. The inner steady state is determined. Around the interior steady state, the deterministic model’s local and global stability is investigated. Furthermore, a stochastic stability study is performed in the neighbourhood of a positive steady state, using analytical estimates of population mean square fluctuations to investigate the system’s dynamics in the presence of Gaussian white noise.

In this work, we looked at a two-zone aquatic habitat where both prey and predators can access the zones. The prey alternates between two zones at random. The growth of prey in the absence of a predator is believed to be logistic in each zone. The inner steady state is determined. Around the interior steady state, the deterministic model's local and global stability is investigated. Furthermore, a stochastic stability study is performed in the neighbourhood of a positive steady state, using analytical estimates of population mean square fluctuations to investigate the system's dynamics in the presence of Gaussian white noise.
In recent years, bio mathematics has been applied to a variety of problems in ecology and epidemiology 1,2 . Policymakers and scientists working in the fields of marine fisheries and other fields debated the economic and social benefits of marine protected zones 3 . Because it is expected that adult or juvenile migration will replace depleted fishing grounds beyond the boundaries of the marine reserve, it may be used as a defensive strategy. Extensive and uncontrolled fishing of marine fish can result in the extinction of various fish species. The development of marine reserves where fishing and other extractive activities are restricted could be one answer to these challenges. Marine reserves not only conserve species within the reserve, but they can help boost fish populations in surrounding areas. Managers must address the interaction between fish and their associated habitats in order for marine reserves to be effective.
Marine reserves are a promising alternative to traditional fisheries management methods such as catch limitations, gear restrictions, and/or effort restrictions. It is supposed to benefit fisheries by conserving spawning populations, providing a safe haven for pre-recruiting fish, and exporting biomass to nearby fishing regions. Nonfisheries benefit of marine reserves include biodiversity and ecological structure protection, biological reference regions, and non-consumption recreational opportunities. The advantages of establishing marine reserves can extend beyond the protection of a single overfished fish stock. Marine reserves can protect the marine environment from damage caused by destructive fishing methods, as well as gives an essential opportunity to learn about marine ecosystems. It aids in the study of species dynamics and the development of management strategies for the conservation of all aspects of a marine community. As a tool for conservation and marine environmental management, marine reserves have yielded numerous benefits 4,5 and has done substantial research on the best management of renewable resources as a fishery, which has a direct link to sustainable development. A dynamic model for fishery resources and showed that the use of diesel-powered trawling may lead to the extinction of predator and prey species if trawling efficiency in the catch of prey species is improved, discussed the economic and biological management of renewable resources [6][7][8][9][10][11] . It is looked at the issue of two rival fish species being harvested at the same time 12 . A model for studying the taxation-based regulation of single-species fisheries was proposed 13 . A method for determining the best harvesting policy for two interdependent populations in a Lotka-Volterra environment was presented 14 . A resource-based competitive system in three species and came up with conditions for the system's persistence and global stability 15 . The uniqueness of limit cycles in a harvested predator-prey system with Holling type III functional response was investigated 16 .
The exploitation of a single species and demonstrated that the time-dependent logistic equation with the periodic coefficient has a single positive periodic solution that is globally asymptotically stable for positive solutions 17 and in a prey-predator fishery, it has been provided a mathematical model of nonselective harvesting 18  www.nature.com/scientificreports/ went on to define taxation as a control device for regulating a prey-predator fishery in their subsequent work 19 . It was also proposed an ideal for studying selective harvesting in a prey-predator fishery, which included a time delay in the harvesting period 20 . A dynamic model created for a single-species fishery that is largely reliant on a supply that is growing logistically 21 . They discovered that when the biomass density of the resource grows, so does the equilibrium density and maximum sustainable output of the fish population. Few researchers developed a harvested population model with diffusional migration recently 22,23 and proposed a mathematical methodology for studying the impact of a reserved zone on an aqua ecosystem's dynamics 24 . The tools for simulating species movement between zones was presented and it has been found that there is an arrangement in which species interact successfully through adaptations to the point that predators become extinct in each patch if adaptations are not made 25 Several mathematical tools to analyze the potential impact of establishing marine sheltered zones on aqua ecosystems, concluding that establishing a marine sheltered zone can result in significant species decline 26 . Many of the benefits associated with marine protected zones have been well investigated, and the field is a hotbed of theoretical ecology and mathematical biology study 27 . The work was presented on the stability of three and four species with migration, bionomic equilibrium, and the best harvesting strategy 28 .
On the basis of the density effect of fish populations, a unique model of fishing effort was created 29 . To characterise a number of widely used fisheries management strategies, they developed new differential equations. The results show that a control parameter (the amount of the effect of fish population size on fishing effort function) controls not only the rate of population equilibrium but also the equilibrium values. They developed a new fishing effort model by contrasting several strategies. In the case of proportional harvesting, increases in β (the noise amplitude) will cause the population to stabilise more quickly. In proportional threshold harvesting, an increase in β will accelerate the convergence to a stable equilibrium point above the equilibrium solution. Due to the periodic nature of proportional harvesting rate, which alters the value of stable equilibrium point, seasonal harvesting strategy has a higher stable equilibrium point than proportional harvesting strategy.
In order for future generations to benefit from fishing, its long-term viability must be ensured through fisheries management. The best harvesting strategy, however, typically maximises a harvester's economically significant objective function, which could result in the extinction of the resource population. As a result, achieving sustainability has proven to be much more challenging than initially anticipated; overfishing is causing catches to decline and fish populations to become increasingly limited. During his research, he developed an efficient harvesting method and formula to maximise the net profit from harvesting.
Constant harvesting and periodic harvesting are the two logistic growth models that have been used. When overcrowding and resource rivalry are taken into account 30 , the logistic growth model is adequate for animal population growth. They wanted to determine which fish harvesting techniques would result in the highest ongoing production. Second, the analysis forecast the ideal harvesting quantity that can guarantee a steady supply of tilapia fish. The optimal harvesting technique for the chosen fish farm, according to findings that were compared between the two strategies, is periodic harvesting. An ongoing seasonal harvesting approach can increase output and accelerate investment payback. Fish farming does not have enough time to rebuild the fish population because of the continual harvesting.
The management of fisheries involves taking into account how harvesting may affect the environment 31 . While fishermen strive to feed an expanding human population, some fish species have been dreadfully dwindling as a result. It's crucial to strike a balance between ecological and economic needs. There were several deterministic models of fishing populations examined. With a constant harvest rate as well as time-dependent harvesting, three different logistic models-a straightforward one, a skewed one with a quadratic term, and one that illustrates the Allee effect-have all been taken into account. The optimal harvest rate for each population density scenario was identified using optimization and numerical computations.
A mathematical model was developed to study the dynamics of a fishery resource system in an aquatic environment that consists of two zones: a free fishing zone and a reserve zone where fishing is strictly prohibited 21 . Biological and bionomic equilibria of the system are obtained, and criteria for local stability, instability and global stability of the system were derived. It was shown that even if fishery was exploited continuously in the unreserved zone, fish populations can be maintained at an appropriate equilibrium level in the habitat. An optimal harvesting policy is also discussed using the Pontryagin's Maximum Principle. Under continuous harvesting of fish species outside the reserved zone, fish population may be maintained at an appropriate equilibrium level. Optimal harvesting policy has been discussed and concluded that high interest will cause high inflation rate. Analyzed a mathematical model to study the dynamics of a fishery resource system with stage structure in an aquatic environment that consists of two zones namely unreserved zone (fishing permitted) and reserved zone (fishing is strictly prohibited) 32 . In this model they introduce a stage structure in which predators are split into two kinds as immature predators and mature predators. It is assumed that immature predators cannot catch the prey and their foods are given by their parents (mature predators). It is also assumed that the fishing of immature predators prohibited in the unreserved zone and predator species are not allowed to enter inside the reserved zone. The local and global stability analysis has been specified. Biological and bionomical equilibrium points of the system were derived. Mathematical formulation of the optimal harvesting policy was given and its solution was derived in the equilibrium case by using Pontryagin's maximum principle. The vital role of reserved zone in aquatic environment for protection of fishery resources from its overexploitation is discussed by several researchers. The use of marine protect area on both biological and economical aspect was studied 33,34 . They investigated the impacts of the creation of Marine Protected Areas (MPAs), in both economic and biological perspectives. The economic indicator was assumed to be optimally managed. The biological indicator was taken as the stock density of the source. The basic fishery model was serving as the benchmark in comparing results with those that are derived from a model of two patchy populations. A crucial characteristic is the migration coefficient which describes biological linkages between protected and unprotected areas. Both economic and biological criteria are enhanced, after introducing a marine protected area, was presented. A prey-predator fishery www.nature.com/scientificreports/ model with the influence of prey reserve were developed and also noted that the fish population maintains at an equilibrium level in absence or presence of predator provided the population in the unreserved area lies in a certain interval 35 . The result obtained was a high interest rate will cause a high inflation rate. Bio-economic model of a prey-predator fishery with protected area and also discussed the system numerically and observed that marine protect area can be used as an effective management tool to improve resource rent under a number of circumstances 36,37 . Bionomic equilibrium and optimal harvesting 39 is one of the interesting and informative studies on fishery models. Mathematical models (including stochastic and diffusive models) are powerful tool to interpret and catch up the dynamics exhibited by fishery model in multi zone environment. One of such model, which is studied by Srinivas et.al 39 is the basic motive for the current work, which we analyzed the dynamics of stochastic system both analytically and graphically and parametric variations with the help of graphical solutions. One of the interesting objectives of the work is the stochastic stability of a prey-predator system in a two-zone aqua environment with logistic growth patterns. It has been prompted us to use analytical tools to investigate the stochastic stability of an aqua environment 29,38 .

Mathematical model
The following system of four nonlinear ordinary differential equations is the stochastic model of a prey-predator fishery resource: The following qualities are presumptions made by the model: It is a prey-predator system in a two-patch habitat. Both regions are accessible to predators and prey. Each patch must have a consistent appearance. The prey is believed to randomly move between the two regions. The growth of prey in each patch is thought to be logistic in the absence of predators. With inherent growth rates α 1 and α 2 carrying capacities related to population size, the predator consumes the food in the area and grows exponentially as the population grows. In patch-1, x(t) , z(t) , E 1 , K , γ 1 , m 1 , q 1 , r , α 1 represents the biomass density of prey species, predator species, effort applied to harvest the fish population, carrying capacity of prey species, equilibrium ratio of prey to predator biomass, prey mortality rate due to predation, catch ability coefficient, intrinsic growth rate of prey species, and intrinsic growth rates of predators, in that order. In patch-2, the symbols y(t) , w(t) , E 2 , L , γ 2 , m 2 , q 2 , s , α 2 for biomass density of prey species, biomass density of predator species, effort put forth to deplete the fish population, carrying capacity of prey species, equilibrium ratio of prey to predator biomass, prey mortality rate due to predation, catchability coefficient, intrinsic growth rate of prey species, intrinsic growth rate of predator species, and catchability coefficient, respectively, are represented. σ 1 , σ 2 represents the migration rates from patch -1 to patch-2 and vice-versa. β i , i = 1, 2, 3, 4 represents the amplitude of noise on the species prey1, prey2, predator1, predator2 respectively.
where δ i j is the Kronecker delta function and δ is the Dirac delta function. The proposed model is a multi-zonal (patch-1, patch-2) ecological environment, where prey and predator populace exist in two patches named as patch-1 and patch-2. There is an interaction among the species of both the patches is migration. Migration of prey species is allowed from patch-1 to patch-2 and vice-versa. The pictorial representation of the proposed model including migration coefficients m 1 and m 2 is presented as a schematic diagram, which is labeled as Fig. 1 is as follows. Similarly, all the attributes involved in the proposed system (1)-4) are described and framed with best fit of values given in Table 1.
Equilibrium analysis without noise. The interior steady state of the system in the absence of noise (i.e.
, where x * , y * , z * , w * are positive solutions of the following equations: After eliminating z * and w * , Eliminating y * from (5) and (6), we get (7) has a unique positive solution if the following inequalities hold.
Also, w * = γ 2 y * , z * = γ 1 x * ; For y * to be positive, we must have Steadiness without noise. Let's pretend that the given system has a single positive equilibrium point P (x * , y * , z * , w * ) , and we're looking at the dynamics of the system around it. In the absence of noise, the characteristic equation of the Jacobian matrix of the system (1)-(4) is given by   (8) has negative real roots. Because the aforementioned conditions are met, the interior equilibrium point P(x * , y * , z * , w * ) is locally stable. Using the Lyapunov theorem, we will now explore the global stability of the interior equilibrium point P(x * , y * , z * , w * ) of the system (1)-(4) in the absence of noise.
To be negative V ′ , we must have A < x < B and, C < y < D which are both provided in the theorem's expression. As a result, if prey populations fall within a given range in the presence of predators and harvesting, they can be maintained at an optimum equilibrium level.
Stability analysis with noise. The effect of random fluctuations in the environment cannot be included in the approximations generated from deterministic models. The stochastic analysis allows us to investigate the dynamics of any natural ecosystem that is subjected to random environmental changes. The parameters of the system bounce about their mean values in the stochastic model (1)(2)(3)(4). As a result, the previously set equilibrium point will now oscillate about the mean state. The random noise is included into the model as additive Gaussian white noise, and thus any system parameter p reduces to p + βξ(t) , where β ∈ R is the noise amplitude and is the ξ(t) Gaussian white noise process. Because the goal of this study is to visualize the dynamics of the system around the interior equilibrium point, we use the perturbation approach to linearize the model.
Then we have, Using the above perturbations, we identify the linear system of the model (1)-(4) as Applying Fourier transform on the linear system, we get the following algebraic system The above system can be represented in the matrix form as Equation (13) can also be written as We now describe some of the basics of the random population function. If the function Y (t) has a zero mean value, then the fluctuation intensity (variance) of its components in the frequency interval [ω, ω + dω] is S Y (ω)dω , where S Y (ω) is spectral density of Y and is defined as If Y has a zero-mean value, the inverse transforms of S Y (ω) is the auto covariance function The corresponding variance of fluctuations in Y (t) is given by and the auto correlation function is the normalized auto covariance as The effect of random environmental fluctuations on specific species can easily be analyzed with either β 1 = 0 or β 2 = 0 or β 3 = 0 or β 4 = 0 . The analytical evaluation of mean square fluctuations is difficult. But the mean square fluctuations can easily be calculated and visualized for different set of parameters. Figures 2, 3, 4, 5, 6 are the stochastic graphs for the proposed model with additive noise on fish population named as Prey-I, Prey-II, Predator-I, Predator-II in two patchy environment. Figure 2 represents the time series evaluation of fish population for the values of noise intensities of 1.5, 1, 1.5 and 1 for all 4 category fish population respectively along with the attributes of Table 1. Figure 2 clearly shows the little (very less) oscillatory behavior exhibited by all fish population in two patches, which says model system undergone less influence at these noise intensities (1.5, 1, 1.5, 1). System influenced by additive noise is notable in this Fig. 2. Figure 3 represents the time series evaluation of populations for the values of noise intensities of 6, 5, 6 and 5 for all 4 category fish population respectively along with the attributes of Table 1. Figure 3 clearly shows good oscillatory behavior exhibited by all fish population in two patches, which says model system undergone great influence at these noise intensities (6, 5, 6 and 5). System influenced by additive noise is notable in this Fig. 3. Figure 4 represents the time series evaluation of populations for the values of noise intensities of 10, 8, 10 and 8 for all 4 category fish population respectively along with the attributes of Table 1. Figure 4 clearly shows more oscillatory behavior exhibited by all fish population in two patches, which says model system undergone great influence at these noise intensities (10, 8, 10 and 8). System influenced by additive noise is remarkable in this Fig. 4. Figure 5 represents the time series evaluation of populations for the values of noise intensities of 80, 70, 80 and 70 for all 4 category fish population respectively along with the attributes of Table 1. Figure 5 clearly shows highly  Table 1. Figure 6 clearly shows more oscillatory behavior exhibited by all fish population in two patches, which says model system undergone great influence at these noise intensities (200, 150, 200 and 150). System influenced by additive noise is remarkable in this Fig. 6.
Delay analysis without noise. We looked at a system with two populations in separate patches. The stability of the interior steady state is investigated, and the system is shown to be stable. The stability results have been shown both analytically and quantitatively. We can also use a delayed model system to account for the prey  Das and Gazi 31,32 have studied delay differential equation models extensively in the study of numerous ecological systems. All of such systems are two-dimensional and three-dimensional. The analytical examination of the system will be difficult to tractable since the current model is a four-dimensional system, and the expression for the delay parameter values will be convoluted for which the system is stable. As a result, we only solve the system numerically.   Figure 7 shows that the time series evaluation dynamics of populations for different initial values 7,10,12,15 . It is clear from the figures that the system exhibits steadiness with in a very short span of time and the system comes to steadiness after certain period of time. They are varying with time for different initial values. Figures 8, 9 shows that the time series evaluation dynamics for both prey populations. And Figs. 10, 11 shows the time series evaluation dynamics of both predator populations for different values of migration rate from parch-1 to patch-2(σ 1 ). It is clear from the figures that the system exhibits steadiness with in a very short span of time and the system comes to steadiness after certain period of time. They may vary with for different initial values. Figures 12, 13 shows that the time series evaluation dynamics for both prey populations. And Figs. 14, 15 shows the time series evaluation dynamics of both predator populations for different values of migration rate from parch-2 to patch-1(σ 2 ). These figures clearly show that the system demonstrates steadiness within a very short length of time and that it attains steadiness after a specific amount of time. They change throughout time depending on the original initial numeric.  Figure 20 depicts the population time series evaluation at a delay value τ 1 = 2.5 . This figure clearly shows that, in contrast to the predator populations of the system, the prey populations of the system exhibit stability within a very short length of time and come to steadiness after a specific amount of time. The beginning values will determine whether they alter over time. Figures 21, 22, and 23 displays phase portrait diagrams for the populations of the proposed system (X(t), Y(t), W(t); X(t), W(t), Z(t); and Y(t), W(t), Z(t)). They might change throughout time depending on the initial settings of numeric. Figure 25 depicts the population time series evaluation at a delay value τ 1 = 12.5 , which clearly shows that, all populations of the system, exhibit oscillatory behavior and unstable nature. Figures 26,27,28,29 presents phase portrait projections plotted among X(t), Y(t), Z(t) and W(t) at a delay value τ 1 = 11.5 , which clearly shows that, limit cycle oscillations of the system and converges to a stable equilibrium point.

Concluding remarks
With prey migration and prey harvesting, the two-patch aquatic habitat for prey and predator is explored in this paper. The growth of the species is considered in a logistical approach. Because the model is four-dimensional, there are 16 mathematical equilibrium points, but the presence of all species and their dynamics is of great ecological significance. As a result, the stability study is limited to the inner steady state. Interior steady state expressions and existence conditions are derived. Routh-Hurwitz criteria are used to determine the system's local stability behaviour in the absence of noise. The Lyapunov theorem is used to investigate the global dynamics of the system in the absence of noise around the interior equilibrium point. The prey population intervals produced show that the system may be kept globally asymptotically stable if both prey populations are within a given range. Additive white noise perturbations following a Gaussian process are used to approximate random environmental disturbances. The higher the population variance, the more unstable the system is, and the lower the population variance, the better the stochastic system's stability behaviour in the respective area. The analytical conclusions may be confirmed using numerical simulations as a future scope of the current analysis. Distinct sets of parameters can display different stability behaviours in local and global aspects, and numerical simulations can easily identify the model's sensitive parameters. The system behaviour with fluctuations in population densities can be simulated and compared for different sets of noise amplitudes. Furthermore, ideal harvesting levels for prey species can be determined in order to promote species expansion and earnings. Delay induced system analysed both analytically and graphically, delay dynamics are identified, analysed and presented well with the help of numerical simulations. www.nature.com/scientificreports/ Figure 20 shows the population time series evaluation at a delay value τ 1 = 2.5 , which clearly shows that, all populations of the system, exhibit stable nature very soon. All Populations-X(t), Y(t), Z(t) and W(t) attains stability very soon and the system is not completely influenced by delay parameter at this value τ 1 = 2.5 . Prey-I, Prey-II shows very little variation at this delay value and Predator-I, Predator-II are not under the influence of delay at this value. Figure 24 depicts the population time series evaluation at a delay value τ 1 = 10.5 , which clearly shows that, all populations of the system, exhibit little oscillatory for some period of time and after certain period of time all the populations X(t), Y(t), Z(t) and W(t) attains stability. Which says that system undergone the good effect of delay, which is exhibited by its oscillatory behavior. Figure 25 depicts the population time series evaluation at a delay value τ 1 = 12.5 , which clearly shows that, all populations of the system, exhibit gradual and consistent oscillatory and unstable nature. At this value of delay τ 1 = 12.5 all the populations X(t), Y(t), Z(t) and W(t) undergone high influence delay and shown their oscillatory behavior. Which says that system undergone the influence of delay highly, which is exhibited by its continued oscillatory behavior. Figures 26, 27  www.nature.com/scientificreports/ the population time series evaluation at a delay value τ 1 = 11.5 , which clearly shows that, limit cycle oscillations of the system and converges to a stable equilibrium point. Figures 2, 3, 4, 5, 6 depicts the time series plots for the stochastic model (1)-(4) System undergone the influence of Gaussian white noise, exhibits its dynamics in the form of chaotic and oscillatory nature in the graphical solutions. Prey-1 and Predator-1 are less oscillatory and Prey-2 and Predator-2 Parameter analysis is carried for influential attributes in the model (1)-(4) and the values of all attributes are from Table 1.
In this paper, we studied the delay dynamics and stochastic dynamics in view of fishery resource management along with migration of prey species in both the zones (patches). This work is also connected to the impact of harvesting and migration (allowed for prey in both patches) on the proposed model. Few examples on prey-1 and prey-2 are Gold band fish/small pelagic fish/small tunas/shrimps. Few examples on predator-1 and predator-2 are shark fish/Dolphin fish/Great white shark/stone fishes etc.